base_dir <- getwd()
files <- list.files(paste0(base_dir,"/data"))
files <- files[grepl("QUANT", files)]
files <- files[grepl("reg", files)]
print(paste0("Found ",length(files)," QUANT files in Dir."))
## [1] "Found 9 QUANT files in Dir."
rois <- sapply(files, function(x) sub("_QUANT\\.tsv", "", x))
print(paste0(length(rois)," ROIs in total"))
## [1] "9 ROIs in total"
get_counts <- function(roi) {
root_path <- base_dir
counts <- fread(paste0(root_path, "/data/for_seurat/", roi, "_qcd_filtered_markers_withids.csv"))
counts$roi <- roi
counts <- counts %>% select(roi, all_of(setdiff(colnames(.), "roi")))
counts
}
get_counts_allmarkers <- function(roi) {
root_path <- base_dir
counts <- fread(paste0(root_path, "/data/for_seurat/allmarkers_", roi, "_qcd_withids.csv"))
counts$roi <- roi
counts <- counts %>% select(roi, all_of(setdiff(colnames(.), "roi")))
counts
}
arcsine_transform <- function(x) {
y <- x / max(x)
asin(sqrt(y))
}
zscorenorm <- function(marker) {
mu = mean(marker)
sd = sd(marker)
(marker - mu) / sd
}
cluster_seurat <- function(codex.obj, res) {
codex.obj <- FindNeighbors(object = codex.obj, dims = c(1:length(markers_selected)), verbose = FALSE)
codex.obj <- FindClusters(object = codex.obj, verbose = FALSE, resolution = res, n.start = 1)
codex.obj
}
get_normalized_counts <- function(roi) {
root_path <- base_dir
counts <- fread(paste0(root_path, "/data/for_seurat/", roi, "_qcd_filtered_markers_withids.csv"))
counts$roi <- roi
counts <- counts %>% select(roi, all_of(setdiff(colnames(.), "roi")))
for_hist_before <- counts %>% select(contains("Median"))
png(paste0("plots/hist_markers_beforenorm_", roi, ".png"), width = 8, height = 8, units = "in", res = 300)
hist.data.frame(for_hist_before)
dev.off()
cols <- grep("Median", colnames(counts), value = T)
for (j in cols) set(counts, j = j, value = arcsine_transform(counts[[j]]))
for (j in cols) set(counts, j = j, value = zscorenorm(counts[[j]]))
for_hist_after <- counts %>% select(contains("Median"))
png(paste0("plots/hist_markers_afternorm_", roi, ".png"), width = 8, height = 8, units = "in", res = 300)
hist.data.frame(for_hist_after)
dev.off()
counts
}
get_normalized_counts_allmarkers <- function(roi) {
root_path <- base_dir
counts <- fread(paste0(root_path, "/data/for_seurat/allmarkers_", roi, "_qcd_withids.csv"))
counts$roi <- roi
counts <- counts %>% select(roi, all_of(setdiff(colnames(.), "roi")))
cols <- grep("Median", colnames(counts), value = T)
for (j in cols) set(counts, j = j, value = arcsine_transform(counts[[j]]))
for (j in cols) set(counts, j = j, value = zscorenorm(counts[[j]]))
counts
}
# Override the default Seurat read function to use median instead of mean
ReadAkoyaCustom_median <- function (filename, type = c("inform", "processor", "qupath"),
filter = "DAPI|Blank|Empty", inform.quant = c("mean", "total",
"min", "max", "std"))
{
if (!requireNamespace("data.table", quietly = TRUE)) {
stop("Please install 'data.table' for this function")
}
if (!file.exists(filename)) {
stop(paste("Can't file file:", filename))
}
type <- tolower(x = type[1L])
type <- match.arg(arg = type)
ratio <- getOption(x = "Seurat.input.sparse_ratio", default = 0.4)
p <- progressor()
p(message = "Preloading Akoya matrix", class = "sticky",
amount = 0)
sep <- switch(EXPR = type, inform = "\t", ",")
mtx <- data.table::fread(file = filename, sep = sep, data.table = FALSE,
verbose = FALSE)
p(message = paste0("Parsing matrix in '", type, "' format"),
class = "sticky", amount = 0)
outs <- switch(EXPR = type, processor = {
p(message = "Creating centroids coordinates", class = "sticky",
amount = 0)
centroids <- data.frame(x = mtx[["x:x"]], y = mtx[["y:y"]],
cell = as.character(x = mtx[["cell_id:cell_id"]]),
stringsAsFactors = FALSE)
rownames(x = mtx) <- as.character(x = mtx[["cell_id:cell_id"]])
p(message = "Creating meta data", class = "sticky", amount = 0)
md <- mtx[, !grepl(pattern = "^cyc", x = colnames(x = mtx)),
drop = FALSE]
colnames(x = md) <- vapply(X = strsplit(x = colnames(x = md),
split = ":"), FUN = "[[", FUN.VALUE = character(length = 1L),
2L)
p(message = "Creating expression matrix", class = "sticky",
amount = 0)
mtx <- mtx[, grepl(pattern = "^cyc", x = colnames(x = mtx)),
drop = FALSE]
colnames(x = mtx) <- vapply(X = strsplit(x = colnames(x = mtx),
split = ":"), FUN = "[[", FUN.VALUE = character(length = 1L),
2L)
if (!is.na(x = filter)) {
p(message = paste0("Filtering features with pattern '",
filter, "'"), class = "sticky", amount = 0)
mtx <- mtx[, !grepl(pattern = filter, x = colnames(x = mtx)),
drop = FALSE]
}
mtx <- t(x = mtx)
if ((sum(mtx == 0)/length(x = mtx)) > ratio) {
p(message = "Converting expression to sparse matrix",
class = "sticky", amount = 0)
mtx <- as.sparse(x = mtx)
}
list(matrix = mtx, centroids = centroids, metadata = md)
}, inform = {
inform.quant <- tolower(x = inform.quant[1L])
inform.quant <- match.arg(arg = inform.quant)
expr.key <- c(mean = "Mean", total = "Total", min = "Min",
max = "Max", std = "Std Dev")[inform.quant]
expr.pattern <- "\\(Normalized Counts, Total Weighting\\)"
rownames(x = mtx) <- mtx[["Cell ID"]]
mtx <- mtx[, setdiff(x = colnames(x = mtx), y = "Cell ID"),
drop = FALSE]
p(message = "Creating centroids coordinates", class = "sticky",
amount = 0)
centroids <- data.frame(x = mtx[["Cell X Position"]],
y = mtx[["Cell Y Position"]], cell = rownames(x = mtx),
stringsAsFactors = FALSE)
p(message = "Creating meta data", class = "sticky", amount = 0)
cols <- setdiff(x = grep(pattern = expr.pattern, x = colnames(x = mtx),
value = TRUE, invert = TRUE), y = paste("Cell", c("X",
"Y"), "Position"))
md <- mtx[, cols, drop = FALSE]
exprs <- data.frame(cols = grep(pattern = paste(expr.key,
expr.pattern), x = colnames(x = mtx), value = TRUE))
exprs$feature <- vapply(X = trimws(x = gsub(pattern = paste(expr.key,
expr.pattern), replacement = "", x = exprs$cols)),
FUN = function(x) {
x <- unlist(x = strsplit(x = x, split = " "))
x <- x[length(x = x)]
return(gsub(pattern = "\\(|\\)", replacement = "",
x = x))
}, FUN.VALUE = character(length = 1L))
exprs$class <- tolower(x = vapply(X = strsplit(x = exprs$cols,
split = " "), FUN = "[[", FUN.VALUE = character(length = 1L),
1L))
classes <- unique(x = exprs$class)
outs <- vector(mode = "list", length = length(x = classes) +
2L)
names(x = outs) <- c("matrix", "centroids", "metadata",
setdiff(x = classes, y = "entire"))
outs$centroids <- centroids
outs$metadata <- md
for (i in classes) {
p(message = paste("Creating", switch(EXPR = i, entire = "entire cell",
i), "expression matrix"), class = "sticky", amount = 0)
df <- exprs[exprs$class == i, , drop = FALSE]
expr <- mtx[, df$cols]
colnames(x = expr) <- df$feature
if (!is.na(x = filter)) {
p(message = paste0("Filtering features with pattern '",
filter, "'"), class = "sticky", amount = 0)
expr <- expr[, !grepl(pattern = filter, x = colnames(x = expr)),
drop = FALSE]
}
expr <- t(x = expr)
if ((sum(expr == 0, na.rm = TRUE)/length(x = expr)) >
ratio) {
p(message = paste("Converting", switch(EXPR = i,
entire = "entire cell", i), "expression to sparse matrix"),
class = "sticky", amount = 0)
expr <- as.sparse(x = expr)
}
outs[[switch(EXPR = i, entire = "matrix", i)]] <- expr
}
outs
}, qupath = {
rownames(x = mtx) <- as.character(x = seq_len(length.out = nrow(x = mtx)))
p(message = "Creating centroids coordinates", class = "sticky",
amount = 0)
xpos <- sort(x = grep(pattern = "Centroid X", x = colnames(x = mtx),
value = TRUE), decreasing = TRUE)[1L]
ypos <- sort(x = grep(pattern = "Centroid Y", x = colnames(x = mtx),
value = TRUE), decreasing = TRUE)[1L]
centroids <- data.frame(x = mtx[[xpos]], y = mtx[[ypos]],
cell = rownames(x = mtx), stringsAsFactors = FALSE)
p(message = "Creating meta data", class = "sticky", amount = 0)
cols <- setdiff(x = grep(pattern = "Cell: Median", x = colnames(x = mtx),
ignore.case = TRUE, value = TRUE, invert = TRUE),
y = c(xpos, ypos))
md <- mtx[, cols, drop = FALSE]
p(message = "Creating expression matrix", class = "sticky",
amount = 0)
idx <- which(x = grepl(pattern = "Cell: Median", x = colnames(x = mtx),
ignore.case = TRUE))
mtx <- mtx[, idx, drop = FALSE]
colnames(x = mtx) <- vapply(X = strsplit(x = colnames(x = mtx),
split = ":"), FUN = "[[", FUN.VALUE = character(length = 1L),
1L)
if (!is.na(x = filter)) {
p(message = paste0("Filtering features with pattern '",
filter, "'"), class = "sticky", amount = 0)
mtx <- mtx[, !grepl(pattern = filter, x = colnames(x = mtx)),
drop = FALSE]
}
mtx <- t(x = mtx)
if ((sum(mtx == 0)/length(x = mtx)) > ratio) {
p(message = "Converting expression to sparse matrix",
class = "sticky", amount = 0)
mtx <- as.sparse(x = mtx)
}
list(matrix = mtx, centroids = centroids, metadata = md)
}, stop("Unknown matrix type: ", type))
return(outs)
}
LoadAkoyaCustom_median <- function (filename, type = c("inform", "processor", "qupath"),
fov, assay = "Akoya", ...)
{
data <- ReadAkoyaCustom_median(filename = filename, type = type)
coords <- suppressWarnings(expr = CreateFOV(coords = data$centroids,
type = "centroids", key = "fov", assay = assay))
colnames(x = data$metadata) <- suppressWarnings(expr = make.names(names = colnames(x = data$metadata)))
obj <- CreateSeuratObject(counts = data$matrix, assay = assay,
meta.data = data$metadata)
coords <- subset(x = coords, cells = Cells(x = obj))
suppressWarnings(expr = obj[[fov]] <- coords)
for (i in setdiff(x = names(x = data), y = c("matrix", "centroids",
"metadata"))) {
suppressWarnings(expr = obj[[i]] <- CreateAssayObject(counts = data[[i]]))
}
return(obj)
}
plot_avg_heatmap <- function(codex.obj, order_manual = FALSE, row_ord = NULL, col_ord = NULL, scale_by = "row") {
mat <- codex.obj@assays$Akoya@scale.data
mat <- t(mat) %>% as.data.frame()
mat$cluster <- codex.obj$seurat_clusters
# Get means for each cluster
smSub <- mat %>%
group_by(cluster) %>%
summarise_all(mean, na.rm = TRUE) %>%
mutate_all(funs(replace(., is.na(.), 0))) %>%
ungroup()
# Get number of cells per cluster for annotation
annoBarValues <- as.data.frame(table(mat$cluster))$Freq
# Create matrix to be used in the heatmap
if (scale_by == "row") {
mat2 <- smSub %>%
select(-c(cluster)) %>% replace(is.na(.), 0) %>%
as.matrix() %>% t() %>% pheatmap:::scale_rows()
} else if (scale_by == "column") {
mat2 <- smSub %>%
select(-c(cluster)) %>% replace(is.na(.), 0) %>%
as.matrix() %>% pheatmap:::scale_rows() %>% t()
} else print("Select a valid argument for scale_by: row or column.")
## Annotation for cluster
ha = HeatmapAnnotation(Cluster = smSub$cluster,
ClusterID = anno_text(smSub$cluster, gp = gpar(fontsize = 12)))
# Create barplot annotation for cluster size for bottom of heatmap
ba = HeatmapAnnotation(CellCount = anno_barplot(annoBarValues,height=unit(2, "cm")))
mat2[is.nan(mat2)] <- 0
colnames(mat2) <- smSub$cluster
if (order_manual) {
if (!is.null(row_ord)) {
Heatmap(mat2,
row_names_gp = gpar(fontsize = 13),
top_annotation = ha,
bottom_annotation = ba,
column_order = col_ord,
row_order = row_ord,
border = TRUE)
} else {
Heatmap(mat2,
row_names_gp = gpar(fontsize = 13),
top_annotation = ha,
bottom_annotation = ba,
column_order = col_ord,
border = TRUE)
}
} else {
Heatmap(mat2,
row_names_gp = gpar(fontsize = 13),
top_annotation = ha,
bottom_annotation = ba,
border = TRUE)
}
}
all_rois <- ldply(rois, get_counts)
fwrite(all_rois, "data/for_seurat/all_rois_18markers.csv")
all_rois <- fread("data/for_seurat/all_rois_18markers.csv")
all_normalized <- ldply(rois, get_normalized_counts)
all_normalized <- all_normalized %>% select(contains("Median")) %>% as.matrix()
colnames(all_normalized) <- sub(": Cell: Median", "", colnames(all_normalized))
all_normalized <- t(all_normalized)
all_normalized <- Matrix(all_normalized, sparse = TRUE)
colnames(all_normalized) <- 1:ncol(all_normalized)
codex.obj <- LoadAkoyaCustom_median(filename = paste0(base_dir, "/data/for_seurat/all_rois_18markers.csv"), type = "qupath", fov = "all_rois")
codex.obj@assays$Akoya@data <- all_normalized
codex.obj <- ScaleData(codex.obj)
VariableFeatures(codex.obj) <- rownames(codex.obj) # Since the panel is small, treat all features as variable.
codex.obj <- RunPCA(object = codex.obj, npcs = 18, verbose = FALSE, approx = FALSE)
codex.obj <- RunUMAP(object = codex.obj, dims = c(1:18), verbose = FALSE)
saveRDS(codex.obj, "data/global_18markers_codex_obj.rds")
markers_selected <- codex.obj@assays$Akoya@counts@Dimnames[[1]]
codex.obj <- FindNeighbors(object = codex.obj, dims = c(1:length(markers_selected)), verbose = FALSE)
codex.obj <- FindClusters(object = codex.obj, verbose = FALSE, resolution = seq(0.1,1.9,0.2), n.start = 1)
saveRDS(codex.obj, "data/codex_obj_18markers_withclusters.rds")
codex.obj <- readRDS("data/codex_obj_18markers_withclusters.rds")
markers_selected <- codex.obj@assays$Akoya@counts@Dimnames[[1]]
UMAP by ROI
DimPlot(codex.obj, group.by = "roi", raster = T)
Clustree
## Clustree takes forever to run, 0.3 was the resolution chosen
tree <- clustree(codex.obj, prefix = "Akoya_snn_res.")
plot(tree)
ggsave(paste0("plots/global_18marker_clustree.png"), height = 10, width = 13, units = "in", dpi = 300)
min_clusters <- 6
treedata <- tree$data
tree_summary <- treedata %>% group_by(Akoya_snn_res.) %>% summarise(sc3_sd = min(sc3_stability) + sd(sc3_stability), n_cluster = n()) %>% filter(n_cluster > min_clusters)
clustering_res <- as.numeric(as.character(tree_summary$Akoya_snn_res.[which.max(tree_summary$sc3_sd)]))
if(length(clustering_res) > 1) clustering_res <- clustering_res[1]
clustering_res <- 0.3
final_nclust <- tree_summary$n_cluster[which.max(tree_summary$sc3_sd)]
treedata$`Selected resolution` <- treedata$Akoya_snn_res. == clustering_res
tree_summary <- treedata %>% group_by(Akoya_snn_res.) %>% summarise(sc3_sd = min(sc3_stability) + sd(sc3_stability), n_cluster = n()) %>% select(Resolution = Akoya_snn_res., `# clusters` = n_cluster)
# kable(tree_summary, align = "r")
ggplot(treedata, aes(x = Akoya_snn_res., y = size, fill = `Selected resolution`)) +
geom_boxplot() +
geom_jitter(width = 0.2) +
theme_bw() +
xlab("Resolution") +
ylab("Cells per cluster") +
ggtitle("Distribution of cluster sizes per resolution")
ggsave("plots/global_18marker_cluster_res_boxplot.png", height = 5, width = 8, units = "in", dpi = 300)
Clustering resolution * 1
codex.obj <- cluster_seurat(codex.obj, res = clustering_res)
DimPlot(codex.obj, raster = T)
ROI composition of clusters
roi_cluster_counts <- table(codex.obj$roi, codex.obj$seurat_clusters)
roi_cluster_counts <- data.frame(roi_cluster_counts)
ggplot(roi_cluster_counts, aes(x = Var2, y = Freq, fill = Var1)) +
geom_bar(stat = "identity", position = "fill") +
xlab("Cluster")
colnames(roi_cluster_counts) <- c("roi", "cluster", "n")
roi_cluster_counts$pct_of_cluster <- NA
for (i in seq_along(roi_cluster_counts$cluster)) {
cluster_total <- sum(roi_cluster_counts$n[roi_cluster_counts$cluster == roi_cluster_counts$cluster[i]])
# print(cluster_total)
roi_cluster_counts$pct_of_cluster[i] <- (roi_cluster_counts$n[i]/cluster_total)*100
}
max_cluster_pct <- roi_cluster_counts %>% group_by(cluster) %>% slice_max(pct_of_cluster)
# write.csv(max_cluster_pct, row.names = FALSE, "data/05-global_18marker_max_cluster_pct.csv")
for (i in seq_along(roi_cluster_counts$cluster)) {
roi_total <- sum(roi_cluster_counts$n[roi_cluster_counts$roi == roi_cluster_counts$roi[i]])
# print(roi_total)
roi_cluster_counts$pct_of_roi[i] <- (roi_cluster_counts$n[i]/roi_total)*100
}
max_roi_pct <- roi_cluster_counts %>% group_by(roi) %>% slice_max(pct_of_roi)
# write.csv(max_roi_pct, row.names = FALSE, "data/05-global_18marker_max_roi_pct.csv")
Heatmap scaled by row
ht <- plot_avg_heatmap(codex.obj, scale_by = "row")
ht <- draw(ht)
row_ord <- row_order(ht)
col_ord <- column_order(ht)
Heatmap scaled by column
ht <- plot_avg_heatmap(codex.obj, order_manual = TRUE, scale_by = "column", row_ord = row_ord, col_ord = col_ord)
ht <- draw(ht)
Heatmap with all markers
Scaled by row
all_rois_allmarkers <- ldply(rois, get_counts_allmarkers) %>% select(!contains("Background"))
fwrite(all_rois_allmarkers, "data/for_seurat/all_rois_allmarkers.csv")
all_normalized_allmarkers <- ldply(rois, get_normalized_counts_allmarkers)
all_normalized_allmarkers <- all_normalized_allmarkers %>% select(contains("Median") & !contains("Background")) %>% as.matrix()
colnames(all_normalized_allmarkers) <- sub(": Cell: Median", "", colnames(all_normalized_allmarkers))
all_normalized_allmarkers <- t(all_normalized_allmarkers)
all_normalized_allmarkers <- Matrix(all_normalized_allmarkers, sparse = TRUE)
colnames(all_normalized_allmarkers) <- 1:ncol(all_normalized_allmarkers)
codex.obj_allmarkers <- LoadAkoyaCustom_median(filename = paste0(base_dir, "/data/for_seurat/all_rois_allmarkers.csv"), type = "qupath", fov = "all_rois")
if (all(codex.obj$cellid == codex.obj_allmarkers$cellid)) {
codex.obj_allmarkers$seurat_clusters <- codex.obj$seurat_clusters
Idents(codex.obj_allmarkers) <- codex.obj_allmarkers$seurat_clusters
}
codex.obj_allmarkers <- suppressMessages(NormalizeData(object = codex.obj_allmarkers, normalization.method = "CLR", margin = 2))
codex.obj_allmarkers <- suppressMessages(ScaleData(codex.obj_allmarkers))
VariableFeatures(codex.obj_allmarkers) <- rownames(codex.obj_allmarkers)
ht_all_row <- plot_avg_heatmap(codex.obj_allmarkers, order_manual = T, col_ord = col_ord, scale_by = "row")
ht_all_row <- draw(ht_all_row)
row_ord_all <- row_order(ht_all_row)
Scaled by column
ht_all_col <- plot_avg_heatmap(codex.obj_allmarkers, order_manual = T, col_ord = col_ord, row_ord = row_ord_all, scale_by = "column")
ht_all_col <- draw(ht_all_col)
Clustering resolution * 2
codex.obj <- cluster_seurat(codex.obj, res = clustering_res * 2)
DimPlot(codex.obj, raster = T)
ROI composition of clusters
roi_cluster_counts <- table(codex.obj$roi, codex.obj$seurat_clusters)
roi_cluster_counts <- data.frame(roi_cluster_counts)
ggplot(roi_cluster_counts, aes(x = Var2, y = Freq, fill = Var1)) +
geom_bar(stat = "identity", position = "fill") +
xlab("Cluster")
colnames(roi_cluster_counts) <- c("roi", "cluster", "n")
roi_cluster_counts$pct_of_cluster <- NA
for (i in seq_along(roi_cluster_counts$cluster)) {
cluster_total <- sum(roi_cluster_counts$n[roi_cluster_counts$cluster == roi_cluster_counts$cluster[i]])
# print(cluster_total)
roi_cluster_counts$pct_of_cluster[i] <- (roi_cluster_counts$n[i]/cluster_total)*100
}
max_cluster_pct <- roi_cluster_counts %>% group_by(cluster) %>% slice_max(pct_of_cluster)
# write.csv(max_cluster_pct, row.names = FALSE, "data/05-global_18marker_max_cluster_pct.csv")
for (i in seq_along(roi_cluster_counts$cluster)) {
roi_total <- sum(roi_cluster_counts$n[roi_cluster_counts$roi == roi_cluster_counts$roi[i]])
# print(roi_total)
roi_cluster_counts$pct_of_roi[i] <- (roi_cluster_counts$n[i]/roi_total)*100
}
max_roi_pct <- roi_cluster_counts %>% group_by(roi) %>% slice_max(pct_of_roi)
# write.csv(max_roi_pct, row.names = FALSE, "data/05-global_18marker_max_roi_pct.csv")
Heatmap scaled by row
ht <- plot_avg_heatmap(codex.obj, scale_by = "row")
ht <- draw(ht)
row_ord <- row_order(ht)
col_ord <- column_order(ht)
Heatmap scaled by column
ht <- plot_avg_heatmap(codex.obj, order_manual = TRUE, scale_by = "column", row_ord = row_ord, col_ord = col_ord)
ht <- draw(ht)
Heatmaps with all markers
Scaled by row
if (all(codex.obj$cellid == codex.obj_allmarkers$cellid)) {
codex.obj_allmarkers$seurat_clusters <- codex.obj$seurat_clusters
Idents(codex.obj_allmarkers) <- codex.obj_allmarkers$seurat_clusters
}
ht_all_row <- plot_avg_heatmap(codex.obj_allmarkers, order_manual = T, col_ord = col_ord, scale_by = "row")
ht_all_row <- draw(ht_all_row)
row_ord_all <- row_order(ht_all_row)
Scaled by column
ht_all_col <- plot_avg_heatmap(codex.obj_allmarkers, order_manual = T, col_ord = col_ord, row_ord = row_ord_all, scale_by = "column")
ht_all_col <- draw(ht_all_col)
Clustering resolution * 5
codex.obj <- cluster_seurat(codex.obj, res = clustering_res * 5)
DimPlot(codex.obj, raster = T)
ROI composition of clusters
roi_cluster_counts <- table(codex.obj$roi, codex.obj$seurat_clusters)
roi_cluster_counts <- data.frame(roi_cluster_counts)
ggplot(roi_cluster_counts, aes(x = Var2, y = Freq, fill = Var1)) +
geom_bar(stat = "identity", position = "fill") +
xlab("Cluster")
colnames(roi_cluster_counts) <- c("roi", "cluster", "n")
roi_cluster_counts$pct_of_cluster <- NA
for (i in seq_along(roi_cluster_counts$cluster)) {
cluster_total <- sum(roi_cluster_counts$n[roi_cluster_counts$cluster == roi_cluster_counts$cluster[i]])
# print(cluster_total)
roi_cluster_counts$pct_of_cluster[i] <- (roi_cluster_counts$n[i]/cluster_total)*100
}
max_cluster_pct <- roi_cluster_counts %>% group_by(cluster) %>% slice_max(pct_of_cluster)
# write.csv(max_cluster_pct, row.names = FALSE, "data/05-global_18marker_max_cluster_pct.csv")
for (i in seq_along(roi_cluster_counts$cluster)) {
roi_total <- sum(roi_cluster_counts$n[roi_cluster_counts$roi == roi_cluster_counts$roi[i]])
# print(roi_total)
roi_cluster_counts$pct_of_roi[i] <- (roi_cluster_counts$n[i]/roi_total)*100
}
max_roi_pct <- roi_cluster_counts %>% group_by(roi) %>% slice_max(pct_of_roi)
# write.csv(max_roi_pct, row.names = FALSE, "data/05-global_18marker_max_roi_pct.csv")
Heatmap scaled by row
ht <- plot_avg_heatmap(codex.obj, scale_by = "row")
ht <- draw(ht)
row_ord <- row_order(ht)
col_ord <- column_order(ht)
Heatmap scaled by column
ht <- plot_avg_heatmap(codex.obj, order_manual = TRUE, scale_by = "column", row_ord = row_ord, col_ord = col_ord)
ht <- draw(ht)
Heatmaps with all markers
Scaled by row
if (all(codex.obj$cellid == codex.obj_allmarkers$cellid)) {
codex.obj_allmarkers$seurat_clusters <- codex.obj$seurat_clusters
Idents(codex.obj_allmarkers) <- codex.obj_allmarkers$seurat_clusters
}
ht_all_row <- plot_avg_heatmap(codex.obj_allmarkers, order_manual = T, col_ord = col_ord, scale_by = "row")
ht_all_row <- draw(ht_all_row)
row_ord_all <- row_order(ht_all_row)
Scaled by column
ht_all_col <- plot_avg_heatmap(codex.obj_allmarkers, order_manual = T, col_ord = col_ord, row_ord = row_ord_all, scale_by = "column")
ht_all_col <- draw(ht_all_col)
